home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / testmatrix / matsignt.r < prev    next >
Text File  |  1994-12-20  |  1KB  |  64 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Matrix sign function of a triangular matrix.
  4.  
  5. // Syntax:    S = matsign ( T ) 
  6.  
  7. // Description:
  8.  
  9. //    Computes the matrix sign function S of the upper triangular
  10. //    matrix T using a recurrence.  
  11.  
  12. //      Adapted from FUNM.  Called by SIGNM.
  13.  
  14. //    This file is a translation of matsignt.m from version 2.0 of
  15. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  16. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  17.  
  18. //-------------------------------------------------------------------//
  19.  
  20. matsignt = function ( T )
  21. {
  22.   local (T)
  23.  
  24.   if (norm(tril(T,-1),"1"))
  25.   {
  26.     error("Matrix must be upper triangular.");
  27.   }
  28.  
  29.   n = max(size(T));
  30.  
  31.   S = diag( sign( diag(real(T)) ) );
  32.   tol = 0;
  33.   for (p in 1:n-1)
  34.   {
  35.     for (i in 1:n-p)
  36.     {
  37.       j = i+p;
  38.       d = T[j;j] - T[i;i];
  39.  
  40.       if (S[i;i] != -S[j;j])    # Solve via S^2 = I if we can.
  41.       {
  42.         # Get S(i,j) from S^2 = I.
  43.         k = i+1:j-1;
  44.         RHS = 0;
  45.         if (all(k)) { RHS = RHS - S[i;k]*S[k;j]; }
  46.         S[i;j] = RHS  / (S[i;i]+S[j;j]);
  47.  
  48.       else
  49.  
  50.         # Get S(i,j) from S*T = T*S.
  51.         s = T[i;j]*(S[j;j]-S[i;i]);
  52.         if (p > 1)
  53.         {
  54.           k = i+1:j-1;
  55.           s = s + T[i;k]*S[k;j] - S[i;k]*T[k;j];
  56.         }
  57.         S[i;j] = s/d;
  58.  
  59.       }
  60.     }
  61.   }
  62.   return S;
  63. };
  64.